function savings_cv_quadrant = generate_results_cross_validation_quadrant(fold,split,reweight,sample,sample_training,sample_testing,sample_wave,target_wave,tcost,covariate,winter_prevuse,SMC,savedir,fbase)

% This function estimates EWM quadrant rules on the training set and evaluates
% cost savings or energy reduction on the testing set. 

% Inputs: 
% (1) fold: permutation id
% (2) split: share of observations in the testing set
% (3) reweight: use density ratio to reweight waves or not
% (4) sample: cross-sectional dataset for the pooled sample
% (5) sample_training: training sample
% (6) sample_testing: testing sample
% (7) sample_wave: set to 0 to indicate full estimation sample
% (8) target_wave: set to 0 to indicate full estimation sample 
% (4) tcost: private marginal cost of implementing the program per household
%  >0 represents cost savings; =0 represents kWh reduction 
% (5) covariate: indicates covariate for analysis: income, size, vintage,
% minimum of baseline consumption, maximum of baseline consumption, or
% standard deviation of consumption
% (6) winter_prevuse: indicates whether baseline consumption is calculated 
% as the mean of consumption in winter months (Jan and Feb) or as the mean
% of specified pre-treatment periods.
% (7) SMC: =1 use social marginal cost; =0 use retail electricity price
% (8) savedir: directory to save the savings table 
% (9) fbase: string used to indicate the propensity score and baseline
% months specification for output table 

% Output: 
% (1) savings_cv_quadrant: a table summarizing the savings point estimate,
% share of households to treat, total savings for the entire state, and 95%
% confidence intervals for the point estimate on the test set

%% Generate grid data and weighted gu for grid search 
% All cross validation quadrant search uses the two-dimensional grid for
% the full sample. This means that the grid is set, however, the number of 
% households within each grid and the gu at each grid depends on
% the data we put into the grid. 

[k1_target_set,k2_target_set,x1_target_set,x2_target_set,grid_x1_target_set,grid_x2_target_set,boo1,boo2,...
    gu_w, gu_target_set,gu_sample_set,dr,dr_hh,Xu_target_set,nw_target_set,n_target_set,Xu_sample_set,nw_sample_set,n_sample_set,x1_sample_set,x2_sample_set,...
    scaled_att_target_set,scaled_att_sample_set,scaled_att_sample_set_weighted,Ind_target_set,Ind_sample_set] ...
    = generate_inputs_cross_validation_quadrant(reweight,sample,sample_training,sample_testing,sample_wave,target_wave,tcost,covariate,SMC);

%% Minimize cost
tic1= tic;

m_W = Inf*ones(k1_target_set,k2_target_set);
v_x1 = grid_x1_target_set;
v_x2 = grid_x2_target_set;

minW = Inf;

for i1=1:k1_target_set
    boo_i1=(x1_target_set - grid_x1_target_set(i1));
     for i2=1:k2_target_set
        boo_i2=(x2_target_set - grid_x2_target_set(i2));
        for sign1 = -1:2:1 %for i = [-1,1]
            for sign2 = -1:2:1
                select = ((sign1.*boo_i1)>=0).*((sign2.*boo_i2)>=0); 
          
                %must meet both threshold conditions to be selected as treated
                W = sum(gu_w.*select);
                %
                if W<m_W(i1,i2)
                   m_W(i1,i2) = W; 
                end
                %
                if (W < minW)
                    min_i1 = i1;
                    min_i2 = i2;
                    min_sign1 = sign1;
                    min_sign2 = sign2;
                    minW = W;        
                end
            end
        end
    end
end

fprintf('Grid Serach takes %.2f sec\n',toc(tic1))

%% treatment status for each household 

in_Ghat = (min_sign1*(x1_target_set - grid_x1_target_set(min_i1))>=0).*(min_sign2*(x2_target_set - grid_x2_target_set(min_i2))>=0);
vhats_n=0;
vhats_p=0;


 %% Export rule parameters for graphs 

switch winter_prevuse % Use month_since_opower=[-21,-1] as pre-treatment months 
    case 1
        s_winter = '_winter';
    case 0
        s_winter = '';
end
if (tcost)>0 % cost_savings
    if (SMC)
        s_cost = 'SMC';
    else
        s_cost = 'PMC';
    end
elseif (tcost)==0
    s_cost='kwh';
end

s_wave=sprintf('%d%d',sample_wave,target_wave);
s_fold=sprintf('%d',fold);
s_split=sprintf('split%.1f',split);

if reweight==0
    filename_coefs=sprintf('coef_cv_quadrant_%s_%s%s_%s_%s_perm%s.mat',covariate,s_cost,s_winter,s_wave,s_split,s_fold);
else
    filename_coefs=sprintf('coef_cv_quadrant_%s_%s%s_%s_%s_fold%s_reweight.mat',covariate,s_cost,s_winter,s_wave,s_split,s_fold);
end
filename = sprintf('%s_%s',fbase,filename_coefs);
fpathf = fullfile(savedir,filename);
fprintf('Saving "%s" to "%s" ... ',filename,savedir); % must use newline symbol "\n" at the end
save(fpathf,'min_i1','min_i2','min_sign1','min_sign2','m_W',...
'in_Ghat','v_x1','v_x2','Xu_target_set','nw_target_set','gu_w','gu_target_set','n_target_set','covariate','vhats_n','vhats_p','minW',...
'dr','dr_hh','boo1','boo2','Xu_sample_set','nw_sample_set','n_sample_set','gu_sample_set','x1_sample_set','x2_sample_set','scaled_att_target_set',...
'scaled_att_sample_set','scaled_att_sample_set_weighted','Ind_target_set','Ind_sample_set','split','s_fold');

fprintf('   Done!\n');

%% Output for table 
percent=sum(nw_target_set.*in_Ghat)/n_target_set;
savings=sum(gu_target_set.*in_Ghat)/n_target_set;
if (tcost)
    total_savings=savings*n_target_set*12;
else
    total_savings=savings*n_target_set*12/1000;
end

savings_cv_quadrant=nan(1,7);
savings_cv_quadrant=array2table(savings_cv_quadrant,'VariableNames',{'rules','covariate','percent','savings','totalsavings','cilb','ciub'}); 

savings_cv_quadrant.rules='quadrant';
savings_cv_quadrant.covariate = {covariate};
savings_cv_quadrant.percent=percent;
savings_cv_quadrant.savings=savings;
savings_cv_quadrant.totalsavings=total_savings;
savings_cv_quadrant.cilb=0;
savings_cv_quadrant.ciub=0;

end